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1.  Introduction 

About  seven  years  ago,  one  of  the  authors  (Puri  [16])  was  confronted  with  the 
following  biological  phenomenon.  At  time  t  =  0,  each  member  of  a  group  of  hosts 
such  as  animals  is  injected  with  a  dose  of  a  specified  virulent  organism  such  as 
viruses  or  bacteria,  which  elicit  a  characteristic  response  from  the  host  during 
the  course  of  time.  This  response  may  be  death,  development  of  a  tumor,  or 
some  other  detectable  symptom.  If  n(t)  denotes  the  number  of  hosts  not  re¬ 
sponding  by  time  t,  the  plot  against  t  of  either  n(t)  itself  or  of  the  proportion 
n(t)/n( 0)  is  known  as  the  time  dependent  response  curve.  These  response  curves 
differ  with  the  dose  and  with  the  type  of  the  organism.  However,  generally 
speaking,  the  larger  the  injected  dose,  the  sooner  the  host  responds.  The  question 
was  raised  as  to  how  one  could  explain  these  observed  response  curves  through  a 
suitable  stochastic  model.  Upon  a  search  of  the  existing  literature  at  the  time, 
it  was  found  that  most  of  the  models  considered  until  then,  were  based  on  the 
hypothesis  of  existence  of  a  fixed  threshold,  namely,  while  the  organisms  are 
undergoing  certain  growth  processes  within  the  host,  as  soon  as  their  number 
touches  a  fixed  threshold  N,  the  host  responds.  In  [16],  this  hypothesis  was 
abandoned;  first,  because  this  hypothesis  is  not  strictly  correct;  second,  it  is  not 
clear  what  value  one  ought  to  choose  for  N  in  a  given  situation,  and  third,  be¬ 
cause  this  hypothesis  makes  the  algebra  unnecessarily  intractable  due  to  the 
involvement  of  the  first  passage  time  problem.  Instead  an  alternative  hypothesis 
originally  suggested  by  Professor  Lucien  LeCam  was  adopted.  Here,  unlike  in  the 
threshold  hypothesis,  the  connection  between  the  number  Z(t)  of  organisms  in  a 
host  at  time  t  and  the  host’s  response  is  indeterministic  in  character.  More  ex¬ 
actly,  it  is  assumed  that  the  value  of  Z(t )  (or  possibly  of  a  random  variable  whose 
distribution  is  dependent  on  the  process  {Z{t)})  determines  not  the  presence  or 
the  absence  of  response,  but  only  the  probability  of  response  of  the  host.  Mathe¬ 
matically,  this  amounts  to  postulating  the  existence  of  a  nonnegative  risk  func~ 
tion  f(x,  t )  such  that 
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(1)  P{host  responds  during  (t,  t  +  r)|not  responded  until  t  and  Z(t)  =  x] 

=  8f(x,  t)r  +  o(r), 

where  5^0  and  /  satisfies  certain  mild  regularity  conditions.  Stochastic  models 
based  on  this  more  appropriate  alternative  hypothesis  have  been  explored  with  a 
reasonable  amount  of  success  in  a  paper  which  appeared  in  the  last  Berkeley 
Symposium  [16]  and  again  in  a  later  paper  connected  with  bacteriophage  repro¬ 
duction  (Puri  [17]).  In  fact,  in  [16]  it  is  assumed  that  the  risk  function  /  de¬ 
pends  not  only  on  Z{t)  but  also  on  the  integral  f*  Z(t)  dr,  which  is  contemplated 

as  a  measure  of  the  amount  of  toxin  produced  by  the  live  bacteria  during  the 
interval  (0,  t)  assuming,  of  course,  that  the  toxin  excretion  rate  is  constant  per 
bacterium  per  unit  time. 

The  above  models  (see  [16],  [17])  apply  to  the  situations  where  the  response 
causing  agents  are  selfreproducing  such  as  viruses,  bacteria,  and  so  on.  A  natu¬ 
ral  question  arises  as  to  how  a  similar  model  based  on  the  alternative  hypothesis 
would  behave  in  situations  where  the  agent  is  not  selfreproducing.  Such  would 
be  the  case  where,  for  instance,  the  agent  is  a  chemical  poison,  insecticide,  or  a 
drug.  This  then  brings  up  the  typical  phenomenon  that  one  faces  in  what  is 
commonly  described  as  quanta l  response  assays.  The  classical  theory  of  quantal 
response  assays  can  be  found  in  books  such  as  Finney  [7],  [8],  Bliss  [6],  and 
others.  One  of  the  purposes  of  this  theory  has  been  to  help  to  arrive  at  an  estimate 
of  the  relative  potency  of  one  drug  against  another  by  using  measures  such  as 
E.D.  50,  the  dose  which  is  just  about  enough  to  cause  response  among,  on  the 
average,  about  50  per  cent  of  the  subjects.  Here,  typically  the  experimenter 
chooses  a  set  of  doses  of  each  drug  and  tests  each  dose  on  a  batch  of  subjects. 
At  the  end  of  the  test,  the  experimenter  records  how  many  of  the  subjects  re¬ 
sponded.  In  order  to  analyze  the  data  so  obtained,  it  has  been  customary  to  make 
the  following  assumptions. 

(i)  For  each  subject  there  exists  a  tolerance  limit  or  a  threshold  level  T.  This 
limit  for  a  subject  is  the  dose  which  will  be  just  sufficient  to  produce  the  re¬ 
sponse,  so  that  the  subject  will  respond  if  z  ^  T  and  will  not  if  z  <  T,  where  z 
is  the  dose  injected. 

(ii)  The  threshold  level  T  is  assumed  to  be  a  random  variable  varying  over  the 
population  of  subjects,  with  a  common  distribution.  Thus,  the  probability  that  a 
randomly  chosen  subject  responds  after  receiving  a  dose  z,  is  given  by 

(2)  P(z)  =  P(T  ^  g ). 

It  is  common  practice  for  the  experimenters  to  use  log  dose  or  x  *=  log  z, 
known  as  the  dose  metameter  of  z.  Now  if  g(y)  is  a  probability  density  function 

so  that  J°°  g(y)  dy  =  1,  the  form  of  the  distribution  of  log  T  typically  can  be 

represented  by  the  density 

(3)  dQ{x)  =  T)g(y  +  yx)  dx, 
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where  y  and  rj  are  the  usual  location  and  scale  parameters,  respectively.  With  this, 
one  easily  obtains 

(4)  P{z)  =  f^vlos*g(y)  dy. 

In  practice,  the  choice  of  g(y )  and,  hence,  of  the  distribution  of  the  tolerance 
limit  T  is  rather  arbitrary.  Some  of  the  typical  choices  of  g(y)  that  have  been 
used  in  literature  are  given  below: 


(5) 

g(y)  =  (2tt)-^  exp  {-y*/ 2}, 

—oo  <  y  <  oo, 

(6) 

,  s  fsin2y, 

0  ^  y  ^  tt/2, 
otherwise, 

(7) 

g(y )  =  \  sech2  y, 

—oo  <  y  <  oo. 

The  last 

one  has  been  used  by  Berkson  in 

his  well-known  work  in  this  area 

(see  [5]). 

While  the  above  classical  theory  has  been  found  useful  and  is  still  being  used, 
there  are,  however,  certain  unattractive  features  in  it  that  make  one  feel  like 
giving  it  another  look.  Some  of  these  are  as  follows.  First  is  the  above  mentioned 
objection  of  assuming  the  existence  of  a  tolerance  limit  or  a  threshold  level  for 
each  subject  even  though  the  random  element  is  introduced  only  through  allow¬ 
ing  this  limit  to  vary  randomly  from  subject  to  subject.  Second,  the  model  as  it 
stands  does  not  lend  itself  to  the  consideration  of  any  biological  mechanism 
going  on  within  the  host  leading  to  its  response.  And  third,  it  does  not  allow  the 
consideration  of  the  time  when  the  response  actually  occurs  if  it  does;  all  it 
considers  is  whether  the  response  does  or  does  not  occur  within  a  fixed  length 
of  time.  These  same  features  also  underlie  the  more  recent  work  of  Ashford  [1], 
Ashford  and  Smith  [2],  and  Plackett  and  Hewlett  [13],  [14],  in  the  case  of  mix¬ 
ture  of  drugs. 

It  is  the  purpose  of  this  paper  to  give  the  classical  theory  of  quantal  response 
assays  a  fresh  look  and  to  construct  new  stochastic  models  that  attempt  to 
eliminate  the  above  objections.  This  has  been  achieved  by  adopting  the  alter¬ 
native  approach  of  the  nonthreshold  type  as  discussed  above.  For  the  biological 
phenomenon  under  consideration,  a  typical  stochastic  model  of  the  present  type 
would  involve  the  consideration  of  the  following  three  main  components. 

The  input  process.  This  describes  the  manner  in  which  the  drug  is  introduced 
into  the  subject.  We  call  it  the  “input  process.”  One  could  visualize,  depending 
upon  the  situation  in  question,  several  possibilities  of  inputs  such  as  a  continuous 
time  deterministic  input,  discrete  time  deterministic  input,  or  a  random  input 
according  to  some  random  mechanism. 

The  release  process.  This  describes  the  manner  in  which  the  subject  attempts 
to  reduce  the  level  of  the  drug  within  its  body.  This  may  be  carried  out  either 
through  the  process  of  direct  elimination  of  the  drug  through  natural  means  or 
by  changing  the  composition  of  the  drug  itself  through  biochemical  processes. 
We  call  this  the  “release  process.”  In  principle,  this  would  involve  the  mechanism 
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going  on  within  the  body  of  the  subject  that  takes  into  account  the  manner  in 
which  the  subject  copes  with  the  drug.  In  experimental  situations,  the  input  proc¬ 
ess  is  generally  controlled  by  the  experimenter.  The  release  process,  on  the  other 
hand,  is  much  more  involved.  It  requires  a  great  deal  of  experience  and  knowledge 
of  the  biological  system  on  the  part  of  the  experimenter.  It  involves,  in  general, 
a  considerable  amount  of  experimentation  while  probing  into  the  nature  of  the 
release  mechanism.  There  has  been,  in  fact,  much  work  done  in  the  past  in  an 
attempt  to  describe  this  mechanism  for  certain  situations.  For  instance,  the 
compartment  models  of,  among  others,  Sapirstein,  Vidt,  Mandel,  and  Hanusek 
[22] ,  Bellman  [3] ,  [4] ,  are  attempts  toward  a  better  understanding  of  functioning 
of  specific  organs  and  of  various  biological  systems.  Unfortunately,  not  too  many 
of  these  models  are  stochastic  in  nature.  Again,  the  models  in  dam  theory  (see 
Moran  [11],  Gani  [9],  and  Prabhu  [15],  to  cite  only  three  references  from  this 
vast  literature)  could  be  found  suitable  for  combining  the  aspects  of  both  the 
input  and  the  release  processes. 

The  risk  junction.  The  most  important  aspect  that  does  not  appear  to  have 
been  considered  before  in  the  context  of  the  classical  quantal  response  assays  is 
the  consideration  of  a  risk  function  that  relates  the  input  and  the  release  processes 
of  the  drug  to  the  causation  of  the  subject’s  response.  Whether,  in  any  given 
situation,  the  risk  function  depends  only  on  the  level  Z(t )  of  the  drug  at  time  t, 
or  on  some  other  factors  characterising  the  biological  mechanism  going  on  within 
the  body  of  the  subject,  would  entail  a  considerable  knowledge  of  the  biological 
system. 

In  the  next  few  sections,  we  shall  attempt  to  incorporate  the  three  aspects 
listed  above  into  a  stochastic  model.  Although,  this  has  been  done  here  under 
rather  simplified  assumptions,  the  results  do  indicate  that  there  is  something 
to  be  gained  by  approaching  this  problem  from  a  structural  point  of  view.  In 
this  context  the  reader  may  also  find,  among  others,  the  work  of  Neyman  and 
Scott  [12]  of  great  interest.  Here  the  response  causing  agent  is  urethane,  while 
the  response  is  the  appearance  of  tumors  in  mice. 

2.  A  stochastic  model  based  on  a  quantal  response  process 

2.1.  Assumptions  and  notation.  As  a  first  attempt,  we  consider  here  a  simple 
stochastic  model  along  the  lines  discussed  above.  More  comprehensive  models 
incorporating  detailed  mechanisms  suitable  for  certain  situations  shall  be  re¬ 
ported  elsewhere.  Following  the  lines  of  classical  quantal  response  assays,  we 
assume  that  for  each  subject  the  experiment  starts  with  the  administration 
(input)  of  a  single  dose  Z(0)  =  z  at  time  t  =  0,  with  no  other  inputs  thereafter. 
Thus,  if  Z(t)  denotes  the  amount  of  drug  present  at  time  t  in  the  body  of  the 
subject,  it  is  evident  that  with  probability  one  Z(t)  is  nonincreasing  with  t.  The 
release  process  is  assumed  to  have  two  components.  The  first  one  determines  how 
often  and  at  what  times  the  releases  occur,  while  the  second  one  associates  with 
each  such  occurrence  a  nonnegative  random  variable  Y  denoting  the  amount 
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of  the  drug  to  be  released  if  available.  More  specifically,  if  N(t )  denotes  the  num¬ 
ber  of  releases  occurring  during  (0,  t\,  we  assume,  for  simplicity,  that  N(t)  is  a 
Poisson  process  with  parameter  p  >  0.  Also,  given  N(t),  let  Fi,  F2,  •  •  •  ,  Fa(/) 
denote  the  random  amounts  to  be  released  if  available,  at  the  release  time  points 
as  determined  by  the  Poisson  process.  In  particular,  it  is  assumed  that  condi¬ 
tionally  given  N(t),  the  random  variables  Yi,  F2,  •  •  •  ,  Ynu)  are  independently 
distributed  with  a  common  distribution  having  the  probability  density  function 


(8) 


h(y)  =  gexp 


y>  o, 

elsewhere, 


where  /3  >  0.  Of  course,  if  at  any  time  the  random  amount  F,  is  greater  than  the 
amount  actually  available,  all  the  available  amount  is  then  released.  From  the 
above  construction,  it  follows  that 


(9)  z(t)  3=  max^O,  Z(0)  —  £  Fyj>  t  ^  0, 

where  by  convention  F0  =  0.  Under  the  Poisson  process  assumption,  it  is  clear 
that  how  often  and  at  what  times  the  releases  occur  is  not  influenced  by  the 
changes  over  time  in  the  amount  of  the  drug  actually  present.  This,  however, 
may  not  be  realistic  in  certain  situations.  In  Section  5,  we  shall  consider  briefly 
a  more  general  model  incorporating  this  dependence  in  an  appropriate  manner. 
Finally,  we  introduce  what  we  shall  call  a  quanta l  response  process  x(t)  defined  as 


(10) 


if  the  subject  does  not  respond  until  t, 
otherwise, 


where  x(0)  =  1.  Also  it  is  assumed  that 

(11)  P{X(t  +  r)  =  0|X(0  =  1,  Z{t)  =  x]  =  8f(x,  t)r  +  o(r), 

where  5  >  0,  and  /( • ,  • )  defined  for  x  ^  0  and  t  ^  0,  is  a  nonnegative  bounded 
function,  assumed  to  be  continuous  almost  everywhere  with  respect  to  both  of 
its  arguments.  Using  a  standard  argument,  it  is  easy  to  show  that 

(12)  P{x(t)  =  1|«}  =  exp  {-5  f*f(Z(r,  o>),  t) 


where  Z(r,  co)  denotes  the  state  of  the  process  {Z(t)}  at  time  r  for  a  given  realiza¬ 
tion  co  of  this  process.  From  (12),  we  easily  obtain  the  transform 


(13)  E[x(t)  exp  {-sZ(t)}]  =  E  [exp  {  — sZ(<)  -  8  J*f(Z(r),  r)  dr}]> 


where  Re(s)  ^  0.  In  particular,  this  yields 

(14)  P{L  >  <}  =  P{x(t)  =  1}  =  E[x(l)]  =  E  [exp  {-8  j‘f(Z(r),  r)  dr}], 


where  L  is  the  length  of  time  the  subject  takes  to  respond.  Taking  8  in  (14)  as  a 
dummy  variable,  it  follows  that  the  response  time  distribution  can  be  studied 

equivalently  by  obtaining  the  distribution  of  the  integral  JQ  f[Z(r),  r]  dr.  The 
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reader  may  find  this  particular  connection  explored  in  detail  elsewhere  (see  Puri 
[18],  [19],  [20],  [21]).  Again,  in  general,  the  random  variable  L  may  not  be  a 
proper  random  variable,  so  that 

(15)  P{L  —  o o}  =  P{no  response} 

=  lim  E[X(t)]  =  E  [exp  {  -  5  fj  f(Z(r),  r)  dr}]- 

At  this  point,  we  introduce  the  following  notation : 

Wi(t,z,x)  =  P{Z{t )  ^  x,  x(t)  =  1\Z(0)  =  2,  x(0)  =  1}, 

W1(t)  z)  =  P{xW  =  1|Z(0)  =  x(0)  =  1}, 

W(t,  z,  x)  =  P{Z(t)  ^  x\Z(0)  =  z}, 

(16)  0i(0,  z,  x)  =  exp  {—6t}Wi(t,  z,  x)  dt, 

0i(0,  z)  =  Jq  exp  {— Qt)W\{t,  z)  dt, 

0(0,  z,  x)  =  f;  exp  {— 9t}W(t,  z,  x)  dt, 
where  Re(0)  >  0,  0  ^  x  ^  z,  and 

,17v  W(t,0,x)  =  W(t,z,z)  =  1,  x^O, 

^  '  W(t,  z,  x)  =  W\(t,  z,x)  —  0,  x  <  0. 

Here  the  last  line  follows  from  the  fact  that  zero  is  an  absorption  state  for  the 
process  Z{t). 

In  the  next  subsection,  we  shall  attempt  to  obtain  expressions  for  the  quanti¬ 
ties  defined  above,  by  setting  up  the  usual  Kolmogorov  backward  integral  equa¬ 
tions  involving  these  quantities. 

2.2.  Certain  integral  equations  and  their  solutions.  Unless  stated  to  the  con¬ 
trary,  we  assume  henceforth  that  the  risk  function  /  does  not  explicitly  depend 
on  time  t  and  that  it  depends  only  on  the  level  Z{t).  Moreover,  it  is  assumed  that 
fix)  is  differentiable  for  all  x  ^  0. 

By  considering  the  moment  of  the  first  release  during  (0,  t)  and  the  amount  to 
be  released,  it  is  easy  to  establish  the  following  Kolmogorov  backward  integral 
equation  for  the  probability  W\{t,  z,  x)  for  x  <  z, 

(18)  Writ,  z,  x) 

=  nP  f0l  exp  {-[u  +  8f(z)]u}  [f*~x  Wi(t  -  u,z-  y,x)  exp  {-py}  dy 
+  fz  Wi{t  -  u,z  -  y)  exp  {-py}  dy 

Jz—x 

+  exp  (—5/(0)  (t  —  u)}  exp  {—Py}  dy ]  du. 
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Taking  the  Laplace  transform  of  both  sides  of  (18),  we  have  for  Re  0  >  0, 

(19)  [m  +  0  +  8f(z)]<f> i(0,  2,  x)  =  nP  J*  X  4>i(B,  z  —  y,x)  exp  {-py}  dy 

+  up  fe  (kid,  z  -  y)  exp  {-py}  dy 

j  z—x 

+  exp  {  Pz}  [m/(0  +  5/(0))]. 

Similarly,  we  have  the  corresponding  equation  for  Wi{t,  z)  given  by 

(20) 

Wi{t,  z)  =  exp  {-[m  +  mm 

+  m/3  J*  exp  {-[m  +  5/(2) ]w)  J‘  Wi(t  —  u,z  —  y)  exp  {-p y}  dy  du 

+  m  exp  {-/ Qz]  f*  exp  {-[m  +  5 f{z)]u  -  5/(0)  (J  -  w)}  du, 
or  equivalently  in  terms  of  its  Laplace  transform,  by 

(21)  [m  +  0  +  5/(2)] (0,  2)  =  1  +  exp  {-/32>[m[0  +  5/(0)]-!] 

+  m/3  exp  {-182}  j*  exp  {Pv}<h(9,  v )  dw. 

Equation  (21)  can  easily  be  converted  into  the  differential  equation 

(22)  0(  +  0i[/S  +  (5 f  -  m/3)(m  +  6  +  5/)-1]  =  Pin  +  6  +  8f)~\ 

where  <f>[  and  /'  are  the  corresponding  derivatives  with  respect  to  2.  Solving  (22) 
subject  to  the  initial  condition 

(23)  0i(0,O)  =  [0+5/(O)]->, 
we  obtain 

(24)  0i (0,  2)  =  [0  +  m  +  5/(2)  ]~x  exp{-  P  J*  A  (u)  du] 

[(A(0))-'  +  f>  j’  exp  ]p  J”  A(u)  du]  dv 

where 

(25)  A(u )  =  [0  +  5/(w)][0  +  m  +  5 u^O. 

Substituting  (24)  in  (19)  and  solving  (19)  in  an  analogous  manner,  we  have  the 
solution  for  (19)  given  by 

(26)  0i (0j  2,  x)  =  [0  +  m  +  5/(2) ]-1  exp  p  f*  A(u)  du] 

[(A(0))  -1  —exp  {p  J*A(u )  du]  +  P  f*  exp  {«/;  A(u)  du]  dv], 

where  x  <  2.  As  a  check,  letting  x  — » 2  in  (26)  and  subtracting  the  result  from 
(24)  we  obtain,  as  expected, 
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(27)  l;  p{Z®  =  x(0  =  1|Z(0)  =  2,  x(0)  =  1}  exp  {-et}  dt 

=  z)  ~  z}  z— )  =  [0  +  m  +  S/(z)]-1. 

The  expressions  for  the  transforms  given  by  (24)  and  (26),  in  principle,  are 
sufficient  for  determining  the  joint  distribution  of  x(0  and  Z(t).  Unfortunately, 
to  carry  out  the  inversion  of  these  transforms  in  this  generality  is  rather  cum¬ 
bersome.  Later  on,  we  shall  carry  out  their  inversion  for  a  special  case.  Again, 
if  /( 0)  >  0,  using  a  Tauberian  argument,  it  follows  from  (24)  that 

(28)  \p(z)  =  P{L  =  oo|Z(0)  =  z}  =  lim0<£i(0,  z)  =  0, 

0->O 

so  that  L  is  a  proper  random  variable.  In  fact,  using  the  relation 

(29)  0i(0,  z)  —  JQ  exp  {—6t}P{L  >  r}  dt  =  0_1(1  —  A1  [exp  {— 6L}]) 
and  (24),  we  have 

(30)  E[ex p  {— 0L}]  =  1  —  [0  -j-  n  +  5/(z)]-1  exp  /3  f*  A(u )  dwj- 

[fl(4(0))->  +  09  j‘  exp  |^/3  f”  A{u)  dwj-  dv J- 

Now  one  could  easily  obtain  moments  of  L  from  (30).  In  particular,  it  follows 
from  (29)  that 

(31)  E(L)  =  lim  (0,  z)  =  <fo(0,  z). 

8 — >0 

2.3.  Probability  of  no  response  for  the  case  with  /( 0)  =  0.  If  /( 0)  >  0,  this 
would  mean  that  the  response  could  be  caused  even  without  the  presence  of  the 
drug.  However,  in  most  of  the  practical  situations  this  appears  unrealistic,  ex¬ 
cept  when  the  response  is  the  death  of  the  subject.  Even  in  the  latter  case,  one 
could  define  response  as  the  death  caused  by  the  drug  and  not  by  other  causes; 
or  as  an  approximation  to  the  actual  situation  one  could  ignore  the  other  causes, 
in  which  case  /( 0)  =0  would  be  a  reasonable  requirement.  A  more  realistic 
model  of  this  latter  situation  would  be  the  one  which  incorporates  other  causes 
besides  the  one  due  to  the  drug,  since,  in  principle,  all  these  causes  simultaneously 
compete  against  each  other  for  the  life  of  the  subject.  However,  at  present  we 
shall  not  venture  into  this  refinement  and  instead  assume  /( 0)  =0  in  what 
follows.  With  this  assumption,  the  random  variable  L  is  no  longer  a  proper  ran¬ 
dom  variable,  since  the  probability  that  the  subject  never  responds  will  be 
positive.  Again,  in  quantal  response  assays,  where  the  actual  response  times  are 
often  not  reported,  one  is  typically  interested  only  in  the  probability  that  the 
subject  never  responds.  This  is  valid  only  as  an  approximation  assuming,  of 
course,  that  the  subject  has  been  under  observation  for  a  sufficient  length  of 
time.  Using  (24)  with/(0)  =  0,  this  probability,  denoted  by  ^(z),  is  now  given  by 

(32)  \J/(z)  =  P  {subject  never  responds  \Z{ 0)  =  z} 

=  lim  6<t>i(d,  z)  =  +  5/(z)]_1  exp  |-/35  ['  f{u)[n  +  5/(w)]-1  dw}‘ 

0->O  C  yo  J 
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Alternatively,  one  could  easily  justify,  either  from  (28)  and  (21)  with  /( 0)  =  0 
or  by  a  direct  probabilistic  argument,  that  \f/(z)  must  satisfy  the  integral  equation 

(33)  \p{z)  exp  {pz}  [M  +  8f(z)]  =  n  +  nP  J*  exp  {/3v}  \p(v)  dv. 

This,  when  solved  subject  to  the  boundary  condition  \f/( 0)  =  1,  again  yields  the 
expression  (32).  In  the  next  section,  we  exhibit  a  comparison  of  the  quantal 
response  model  with  the  present  one  through  the  use  of  expression  (32). 

3.  A  comparison  of  the  present  model  with  the  classical  one 

It  appears  rather  natural  at  this  stage  to  look  for  some  kind  of  direct  com¬ 
parison  between  the  present  theory  and  the  classical  one.  Unfortunately,  there 
does  not  appear  to  be  any  simple  way  of  making  such  a  comparison,  mainly 
because  the  two  theories  are  based  on  entirely  different  points  of  view.  However, 
if  we  insist  on  making  one  for  the  sake  of  amusement,  the  only  way  which  ap¬ 
pears  reasonable  is  to  equate  the  end  result  common  to  both  the  theories.  More 
specifically,  by  equating  the  probability  of  no  response  under  the  classical  theory, 
namely, 

(34)  1  -  P{z]  =  P{T  >  z}  =  f“+nk"0(y)  dy, 

to  the  probability  of  no  response  under  the  present  model,  namely  rj/(z),  we  ask 
what  risk  function  /( • )  of  the  present  model  would  correspond  to  a  given  density 
function  g(y)  used  in  the  classical  theory.  To  this  end,  one  can  easily  solve  (32) 
for  /(z)  in  terms  of  \f/  yielding 

(35) 

5/(2)  =  [1K2)]-1  exp  {-p z}  [m(1  -  <A(z)  exp  {i 3z})  +  nP  f*  exp  (/3v}^(v)  dvj- 

Now  by  replacing  \p(z)  with  1  —  P(z),  one  obtains  the  desired  risk  function  / 
corresponding  to  a  given  density  g  of  the  classical  model.  For  instance,  the  risk 
function  corresponding  to  the  normal  density  (5)  is  given  by 

(36)  5/(2)  =  /x  exp  {-pz}  [1  —  H{y  +  17  log  z)]~x 

£l  —  exp  {/3z}  (1  -  H(y  +  77  log  2)) 

+  P  j*  exp  {/9v}  (1  -  H(y  +  17  log  v ))  dv]  > 

where  H(x)  =  (2t)  -  ^  Jx  exp  {— r2/2}  dr.  Similarly,  for  the  density  function 
of  (7)  we  have 

(37)  8f(z)  =  n  exp  {- Pz }  [1  +  exp  {2(7  +  77  log  2)}] 

(l  —  exp  {pz}  [1  +  exp  {2(7  +  77  log  z)}]-1 

+  P  j*  exp  {|Sv}  [1  +  exp  {2(7  +  77  log  v)}]_1  dv)- 
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As  expected,  since  P( 0)  =  0,  we  have  /(0)  =  0  in  the  above  formulas.  Similar 
expressions  can  be  obtained  for  /  that  correspond  to  other  densities  often  used 
in  the  classical  theory.  Unfortunately,  as  is  evident,  all  such  expressions  will 
usually  be  complicated,  so  that  there  appears  to  be  no  rationale  for  choosing 
one  or  the  other  form  of  the  risk  function  in  practice.  In  the  next  section,  we 
consider  the  simplest  form  of  the  risk  function,  namely  the  linear  function  8f(x)  = 
8x ,  which  appears  reasonable  at  least  as  a  first  approximation.  The  results  ob¬ 
tained  by  using  this  simple  risk  function  are  then  applied  to  some  observed  data. 


4.  An  application  of  the  model  to  observed  data 

We  shall  now  restrict  ourselves  to  the  case  of  a  linear  risk  function  with  8f(x)  = 
8x.  For  this,  we  have  from  (26), 

(38)  <f>i(d,  z,x)  -  (0  +  m  +  5z)-1[C(2)]“1  exp  {-# z} 

£(1  +  /I0-1)  —  C(x)  exp  {fix}  -f  j 3  f*  C(v)  exp  {fiv}  dv J, 

where  x  <  z  and 

(39)  C(v )  =  [(0  +  ix)(6  +  m  +  to)-1]*'*,  «  ^  0. 

Also  expression  (24)  now  takes  the  form 

(40)  &(0,  z)  =  [(0  +  n)(6  +  m  +  to)"1] [0(h)]-1  exp  {-fiz} 

[0-1  +  (P(n  +  0)-1)  j*  C(v)  exp  {fiv}  <to]- 

The  transforms  (38)  and  (40)  can  be  easily  inverted  to  produce  the  expressions 
for  Wi(t,  z,  x)  and  W\(t,  z),  respectively;  for  instance,  when  fin/ 8  is  not  an  integer, 

(41)  Wi(t,  z)  =  fo|,  [At  (/f  -  l)  (^)‘  exp 

+  Ak  h+iAW  exp  {  (m  +  te)t} ]» 

where 

(42)  Ak(x)  =  x(x  +  l)(z  +  2)  •  •  •  (x  +  k  —  1),  k  ^  1;  A0(x)  m  1, 

and  for  a  >  0, 

(43)  Ik,a(x)  =  exp  {-ay}  dy,  k  ^  1;  I0,a(x)  =  1. 

It  is  easy  to  verify  that  by  letting  8  tend  to  zero  in  (38),  one  obtains 

(44)  lim  0i(0,  z,  x)  =  0(0,  z,  a) 

=  nW  +  m)]-1  exp  {— 00(0  +  m)-1(2  —  x)},  x  <  z, 
a  result  for  the  process  Z(t)  alone  without  the  consideration  of  the  quantal  re- 
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sponse  process  %(£)•  Again,  for  the  present  case  with /(a:)  =  x ,  we  have  from  (32) 
the  expression  for  the  probability  of  no  response,  given  by 

(5  \— i+/W« 

1  +  -  Z) 

It  is  this  probability  which  is  relevant  to  the  fitting  of  our  model  to  appropriate 
data  on  quantal  response  assays.  Expression  (45)  contains  essentially  two  param¬ 
eters,  since  p  and  8  always  appear  as  p/ 8.  However,  for  fitting  the  above  formula 
to  suitable  data,  it  was  found  convenient  to  introduce  the  reparameterization 

(46)  P  =  -»  X  =  p(/ 3  —  p), 

so  that 

(47)  ip(z)  =  (1  +  pzY'*'  exp  j-  (p  + 

Let  Z0  (E.D.  50)  denote  the  dose  which  will  produce  a  response  with  probability 
one  half.  From  (47),  it  follows  that  Zo  satisfies  the  relation 

(48)  log  ^  =  —  ^P  +  Zo  +  ^  loS  (1  4“  pZ o). 

Formula  (41)  was  fitted  to  the  data  based  on  a  study  of  the  toxicity  of  an  insecti¬ 
cide  known  as  Deguelin.  The  data  are  due  to  Martin  [10]  and  have  also  been 
used  by  Berkson  [5]  in  an  attempt  to  fit  the  classical  model  of  the  quantal 
response  assays.  In  the  study  proper,  concentrations  at  different  dose  levels  Zi 
of  Deguelin  were  prepared  in  an  alcohol  medium.  These  were  then  sprayed  on 
groups  of  respective  sizes  n*  of  the  test  insects  (adult  apterousfe  male)  Aphis 
rumicis.  These  sprayings  were  performed  in  a  carefully  controlled  way  using  a 
special  atomizer.  After  spraying,  without  further  handling  the  insects  were  placed 
in  tubes  with  a  small  amount  of  bean  foliage.  They  were  checked  after  about  20 
hours  for  the  number  r»  of  deaths  in  the  fth  group.  These  data  are  given  in  Table 
I. 


TABLE  I 

Martin’s  Data  on  Toxic  Effect  of  Deguelin 


Concentration 
mg/litre  Zi 

10.1 

20.2 

30.3 

40.4 

50.5 

Total  number  n» 

48 

48 

49 

50 

48 

Number  of  deaths  r< 

18 

34 

47 

47 

48 

Formula  (47)  was  fitted  to  the  above  data  by  using  the  standard  method  of 
minimum  chi  square.  The  fit  appears  quite  satisfactory,  since  the  observed  value 
of  the  chi  square  is  3.69  (three  degrees  of  freedom),  which  is  not  significant  at 
the  five  per  cent  level,  where  the  table  value  is  7.81.  Also,  the  method  of  maxi¬ 
mum  likelihood  led  to  the  estimates  for  the  parameters  p  and  p,  along  with  their 
standard  errors,  as  given  in  (49).  Using  these  and  the  relation  (48),  an  estimate 
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20  of  E.D.  50  was  obtained  by  using  a  computer  search  procedure  for  finding 
the  appropriate  root  of  (48).  The  various  estimates  of  the  standard  errors, 
as  given  here,  are  based  on  the  standard  large  sample  formulas  valid  for  the 
maximum  likelihood  estimates: 

X  =  0.00526,  p  =  0.02428,  =  13.117 

S.E.(X)  =  0.00102,  S.E.(p)  =  0.0143,  S.E.(Z0)  =  1.435. 

As  a  passing  remark,  it  may  be  appropriate  to  mention  here  that  we  also  fitted 
Formula  (47)  to  data  reported  in  [5]  which  pertain  to  responses  to  certain  bac¬ 
teria.  Here,  as  expected,  the  fit  was  considerably  worse.  For  the  four  degrees 
of  freedom  available  in  that  case,  the  observed  chi  square  was  12.9.  This  being 
significant  indicates  the  sensitivity  of  the  present  model  to  situations  where  the 
response  causing  agent  is  selfreproducing.  The  models  appropriate  for  such  situ¬ 
ations  have  already  been  dealt  with  elsewhere  (see  [16],  [17]).  The  present  model 
is,  of  course,  not  designed  for  such  situations. 


5.  A  model  with  a  generalization  of  the  release  process 


In  the  release  process  as  adopted  in  the  above  model,  how  often  and  at  what 
times  the  releases  occur  is  not  influenced  by  the  changes  over  time  in  the  amount 
of  the  drug  actually  present  in  the  subject.  When  the  above  work  was  presented 
at  this  Symposium  a  question  was  raised  from  the  audience  inquiring  whether  it 
was  possible  to  modify  the  release  process  of  the  model  in  order  to  take  into 
account  the  possible  effect  of  the  changes  over  time  in  the  level  of  the  drug  on 
the  frequency  of  the  releases.  We  attempt  here  to  accomplish  this  through  a 
generalization  of  the  Poisson  process  of  Section  2.  Let  n(z)  be  a  nonnegative 
bounded  function,  which  may  be  called  the  risk  function  for  the  release,  such 


that 

P{a  release  occurs  during  (t,  t  +  r)\Z(t)  =  z}  —  h(z)t  +  o(t), 

^  ’  P{more  than  one  release  occurs  during  ( t ,  t  +  r)|Z(i)  =  z}  =  o(t). 


The  random  variables  Yh  F2,  •  •  •  ,  denoting  the  amounts  to  be  released,  if  avail¬ 
able,  at  the  release  points  governed  by  (50),  are  as  before  independently  dis¬ 
tributed  with  the  common  distribution  given  by  (8).  Clearly,  when  ju(z)  is  a 
positive  constant,  we  are  back  to  the  case  of  the  Poisson  release  process.  All  the 
other  assumptions  of  the  model  as  outlined  in  Section  2  remain  the  same  with 
the  exception  of  (50);  we  assume  that/(0)  =  0.  It  may  be  remarked  here  that 
there  is  no  loss  in  generality  so  far  as  the  distribution  of  the  quantal  response 
process  x(0  is  concerned,  if  we  allow  ju(0)  to  be  positive.  In  the  latter  case  we 
can  still  talk  fictitiously  of  the  releases,  even  though  the  level  of  the  drug  may  be 
zero.  Thus,  we  assume  that  ju(0)  >  0,  for  convenience.  Let  N(t)  denote  the  num¬ 
ber  of  releases  occurring  during  (0,  t].  Also  we  introduce  the  following  notation: 


Vi(k,  t,  z)  =  P{x(t)  =  1,  N(t)  =  k\Z(0)  =  2,  x(0)  =  1),  k  ^  0, 
=  P(x(0  =  MZ(0)  =  z,  x(0)  =  1}, 

V(k,  t,  z)  =  P{N(t)  =  k\Z(0)  =  2}, 


(51) 


k  ^  0. 
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It  is  not  too  difficult  to  show  that  the  random  variable  N(t)  is  a  proper  random 
variable  for  every  t  ^  0,  so  that  the  probability  of  an  infinite  number  of  releases 
occurring  during  a  finite  time  interval  is  zero.  As  such 

(52)  Vx{t,z)  =  i  Vi(k,t,z). 

k  =  0 

Again  taking  into  account  the  first  release,  if  it  occurs,  it  is  easy  to  establish  the 
following  system  of  recurrence  relations  for  the  F : 

(53)  Fi(0,  t,  z )  =  exp  {-[n(z)  +  8j(z)]t}, 

(54)  V1(k,t,z) 

=  i t(z)  Jq  exp  {-  [ju (z)  +  Sf(z)]u}  [jj  p  exp  {-py}  Fi (k  -  1,  t  -  u,  z  -  y)  dy 
4-  exp  {— (3z}V(k  —  1,  t  —  u,  0)J>  k  ^  1. 


V* (k,  6,z)  =  f  exp  {— 0f}  Vi  (k,  t,  z )  dt, 

(55)  7„ 

V*(k,  0,  z)  =  I  exp  {— 9t}V(k,  t,  z)  dt, 

where  Re  0  >  0.  Then  from  (53)  and  (54),  we  have 

(56)  FT(0  ,9,z)  =  [9  +  »(z)  +  8f(z)]~\ 

(57)  Vi(k,  9,  z)  =  n(z)[n(z)  +  8f(z)  +  0]_1  [exp  {~Pz}V*(k  -  1,  0,  0) 

+  P  f“  exp  {-/ 8y}Vt(k  -  1,  0,  z  -  y)  dy ]>  k  ^  1. 


Clearly, 
(58) 
so  that 


F(M,0)  =  {-M(0)0, 


(59)  V*(k,  e,  0)  =  D<(0)]‘W0)  +  «]-*-•. 

Using  this,  one  can  easily  solve  the  system  (56)  and  (57)  recursively.  However, 
our  aim  is  to  obtain  \f/(z),  the  probability  of  no  response.  To  this  end,  summing 
(56)  and  (57)  over  the  possible  values  of  k,  we  obtain 

(60)  Vt(9,z) 

=  [0  +  n(z)  +  S/(z)]_1  [l  4-  m(z)  exp  {-0z}  {0~x  4 -  P  f*  exp  {pv)V\(9,  v)  dv}  ]• 

We  assume  now,  for  simplicity,  that  besides  f(z),  the  risk  function  y(z)  is  also 
differentiable  for  z  ^  0.  With  this  (60)  can  be  easily  transformed  into  the  differ¬ 
ential  equation 

dF? 

(61)  4-  Wf  4 -P6  +  p8f)  -  n'(6  4-  5/)][m(0  4-  m  4-  m~lVt 

=  ^  ~  ~  j  [0  4-  m  4-  0/]“S 
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where  f  and  y  denote,  respectively,  the  derivatives  of  /  and  y.  Here  we  have 
suppressed,  for  convenience,  the  arguments  of  all  the  functions  such  as/,  y,  and 
so  forth.  Equation  (61)  can  be  solved  easily  subject  to  the  initial  condition 
Vt(d,  0)  =  1/0,  yielding 


(62)  Vt(d,z) 


y(z){6  +  mo) 
yo[9  +  y(z)  +  5/(z)] 


exp  {-/3 B(z)} 


Mo 

0  +  Mo 


y'(s))(y(s))~2  exp  {(3B(s)}  ds  > 


where  yo  =  y(0)  and 
(63) 


t>(  \  _  [*  a  +  8f(v) 

W  ~  Jo  9  +  y(v)  +  8f(v) 


dv. 


Finally,  since  ${z)  =  lim  9V*(9,  z),  it  follows  from  (62),  that 

0->O 


(64)  $(z)  =  y(z)[y(z)  +  5/(z)]-1  exp  |-/3  J“  8f(v)[y(v )  +  8f(v)]  1  dv }• 


This  then  is  the  generalization  of  formula  (32),  where  y(z)  is  assumed  to  be  a 
positive  constant.  Finally,  for  the  special  case  with  f(x)  =  x  and  y(x)  =  y0  +  vx 
such  that  y0  +  vx  >  0  for  0  ^  x  ^  z,  we  have 


(65)  *{z)  =  exp  .}(l  +  2)(l  +  .) 


—  l+|88/»0(i+v)-* 


2^0. 


6.  Discussion 

The  present  work  is  inspired  by  an  earlier  work  of  one  of  the  authors  (Puri 
[16])  and  by  the  need  of  giving  a  fresh  look  to  the  classical  theory  of  quantal 
response  assays  (see  Finney  [7]),  which,  in  the  opinion  of  the  authors,  appears 
to  have  certain  unappealing  features.  Most  of  the  mathematical  models  of  ran¬ 
dom  phenomena  incorporate  assumptions,  which  tend  to  simplify  the  real  situ¬ 
ation,  yet,  by  now,  it  is  evident  that  there  are  certain  fundamental  differences 
in  the  approach  adopted  here  from  the  one  classically  used.  For  instance,  the 
present  approach  permits  the  consideration  of  the  response  time,  while  the  clas¬ 
sical  one  does  not.  Unlike  the  classical  approach,  the  present  one  is  based  on  a 
nonthreshold  hypothesis  which  appears  more  appealing.  Most  importantly,  how¬ 
ever,  the  present  model  allows  ample  room  for  the  consideration  of  the  mech¬ 
anism  of  the  causation  of  the  response,  while  the  classical  theory  does  not.  The 
mechanism  incorporated  in  the  model  studied  here  may  be  oversimplified  for 
certain  situations.  However,  this,  in  general,  can  easily  be  rectified  by  incorporat¬ 
ing  more  complicated  yet  realistic  mechanisms  into  the  present  theory,  usually, 
of  course,  at  the  cost  of  making  the  algebra  more  involved. 

In  the  present  model  the  only  input  allowed  is  at  the  start  of  the  experiment. 
However,  this  can  easily  be  extended  to  cover  the  general  case,  where  the  input 
pattern  over  time  is  controlled  and  determined  ahead  of  time  by  the  experimenter 
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(see  Neyman  and  Scott  [12]).  Also  situations  such  as  exposure  to  natural  radia¬ 
tion,  or  to  specific  chemicals  as  part  of  certain  occupational  hazards,  involve 
perhaps  a  random  mechanism  for  the  input  process.  Such  models  involving  more 
elaborate  input  and  release  processes  are  reported  elsewhere  (see  Senturia  [23]). 

The  classical  theory  of  quantal  response  assays  has  been  extended  to  the  case 
of  multiple  responses  to  one  or  several  drugs  (see  [1])  or  to  the  case  of  a  single 
response  to  a  mixture  of  drugs  (see  [2]).  It  appears  worthwhile  to  examine  and 
extend  the  present  approach  to  cover  these  cases.  Also,  deeper  models  along  the 
present  lines,  while  incorporating  the  role  of  the  defense  mechanism  utilized  by 
the  subject  in  order  to  cope  with  the  drug,  are  very  much  needed.  This  mech¬ 
anism,  of  course,  may  vary  considerably  from  one  situation  to  another.  In  sev¬ 
eral  situations,  to  gain  knowledge  of  this  mechanism  itself  would  need  a 
considerable  amount  of  further  experimentation. 

Again,  in  many  situations,  it  may  appear  realistic  to  consider  the  risk  function 
/  not  only  dependent  on  the  level  Z(t)  of  the  drug  but  also  on  some  other  relevant 
functionals  of  the  process  Z(t).  (See,  for  instance,  [16],  [17],  and  the  work  done 
at  the  Statistical  Laboratory,  University  of  California,  Berkeley  [24],  [25].) 

In  the  present  model,  a  special  form  (8)  of  the  common  distribution  of  Yh  F2, 
•  •  •  ,  the  amounts  released,  was  assumed.  This  can  be  generalized  to  the  case 
with  an  arbitrary  distribution  function,  say  H(y),  for  the  random  variables  Y. 
One  can  easily  set  up  the  integral  equations  analogous  to  (19)  and  (21)  for  this 
case.  For  instance,  equation  (21)  now  takes  the  form 

(66)  [n  -f*  0  -f-  8f(z)]<f>i(d,  z ) 

=  1  +  *.[1  -  ff(*)][0  +  «/(0)]-'  +  11  f‘  Me,  z  -  y )  dH(y). 

Unfortunately,  however,  the  solution  of  these  equations  becomes  relatively 
cumbersome. 

Finally,  it  is  hoped  that,  in  due  course,  the  approach  adopted  here  will  find 
its  proper  place  in  its  usefulness  in  comparison  to  the  classical  approach.  This 
will  emerge  even  more  when  the  experimenter  wishes  to  use  the  data  on  response 
times  of  the  subjects  for  an  appropriate  analysis,  rather  than  only  on  whether 
the  subject  does  or  does  not  respond  in  a  given  period  of  time. 


7.  Summary 

The  classical  theory  of  quantal  response  assays  (see  [7])  is  based  on  the  hy¬ 
pothesis  of  existence  of  a  threshold  level  T  (tolerance  limit),  such  that  if  the 
injected  dose  z  of  the  drug  is  smaller  than  T,  the  subject  does  not  respond,  and 
it  does  respond  if  z  T.  The  threshold  level  T  is  assumed  to  vary  randomly 
over  the  population  of  subjects  with  a  common  distribution,  usually  with  an 
arbitrarily  chosen  form.  As  it  stands,  the  classical  theory  has  several  unattractive 
features.  First,  the  hypothesis  of  existence  of  a  threshold  level  may  not  be  strictly 
correct.  Second,  the  classical  model  does  not  lend  itself  to  the  consideration  of 
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any  biological  mechanism  leading  to  the  subject’s  response.  And  third,  it  does 
not  allow  for  consideration  of  the  time  when  the  response  actually  occurs,  if  it 
does.  In  view  of  these  objections,  the  present  paper  gives  a  fresh  look  to  the 
problem.  As  a  result,  a  new  stochastic  model  is  constructed  along  the  more 
realistic  lines  (see  [16],  [17])  adopted  elsewhere  for  a  similar  situation.  Here  it  is 
assumed  that 

(67)  P{X(t  +  r)  =  0|x(<)  =  1,  Z{t)  =  x}  =  8f(x)r  +  o(r), 

where  x(0  is  one  if  the  subject  has  not  responded  until  time  t,  and  is  zero  other¬ 
wise;  the  random  variable  Z(t)  denotes  the  level  of  the  drug  at  time  t;  8  >  0  and 
/  is  a  nonnegative  risk  function  defined  for  all  x  ^  0  and  is  assumed  to  satisfy 
certain  mild  regularity  conditions.  The  process  {Z(t)}  is  assumed  to  involve  a 
certain  random  release  mechanism.  Here  Z{ 0)  is  the  dose  administered  at  time 
t  =  0.  Under  these  assumptions,  the  joint  distribution  of  x(0  and  Z(t)  is  studied 
for  an  arbitrary  risk  function  /.  In  particular  with  /(0)  =  0,  the  probability  of 
no  response  is  obtained.  Assuming  f(x)  =  x,  this  probability  is  then  fitted  to  cer¬ 
tain  observed  data.  Finally,  an  indirect  comparison  of  this  model  is  made  with 
the  classical  one. 
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